clc
clear
close all

%%
TimeDuration = 0.1;
SampleRate = 200e3;
MasterClockRate = 600e3;
ModulatorOrder = 4;
SamplePerSymbol = 8;
ModulatorConfig.ModulatorIndex = 0.5;
ModulatorConfig.InitialPhaseOffset = 0;
%
% CarrierFrequency = 200e3;
% IqImbalanceConfig.A = 3;
% IqImbalanceConfig.P = 10;
% PhaseNoiseConfig.Level = -50;
% PhaseNoiseConfig.FrequencyOffset = 20;
% MemoryLessNonlinearityConfig.Method = 'Saleh model';
% MemoryLessNonlinearityConfig.InputScaling = 0;
% MemoryLessNonlinearityConfig.AMAMParameters = [2.1587 1.1517];
% MemoryLessNonlinearityConfig.AMPMParameters =  [4.0033 9.1040];
% MemoryLessNonlinearityConfig.OutputScaling = 0;
% MemoryLessNonlinearityConfig.ReferenceImpedance = 1;
% ThermalNoiseConfig.NoiseTemperature = 290;
% AGCConfig.AveragingLength = 256 * 4;
% AGCConfig.MaxPowerGain = 400;
%
% source = RandomSource(SampleRate=SampleRate, ...
%     TimeDuration=TimeDuration, ...
%     ModulatorOrder=ModulatorOrder, ...
%     SamplePerSymbol=SamplePerSymbol);
%
% x = source();
%
% %% 单天线场景
%
% txRF = TRFSimulator(StartTime=0.2, ...
%     SampleRate = SampleRate, ...
%     CarrierFrequency=CarrierFrequency, ...
%     IqImbalanceConfig=IqImbalanceConfig, ...
%     PhaseNoiseConfig=PhaseNoiseConfig, ...
%     MasterClockRate = MasterClockRate, ...
%     MemoryLessNonlinearityConfig=MemoryLessNonlinearityConfig);
% rayChannel = Rayleigh(CarrierFrequency=CarrierFrequency, ...
%     SampleRate=SampleRate, PathDelays=0, ...
%     AveragePathGains=0, MaximumDopplerShift=0);
%
% ricChannel = Rician(CarrierFrequency=CarrierFrequency, ...
%     SampleRate=SampleRate, PathDelays=0, ...
%     AveragePathGains=0, MaximumDopplerShift=0);
%
% rxRF = RRFSimulator(StartTime=0, TimeDuration=2, SampleRate=SampleRate, ...
%     NumReceiveAntennas=1, CarrierFrequency=200e3, Bandwidth=20e3, ...
%     MasterClockRate=MasterClockRate, ...
%     MemoryLessNonlinearityConfig=MemoryLessNonlinearityConfig, ...
%     ThermalNoiseConfig=ThermalNoiseConfig, ...
%     PhaseNoiseConfig=PhaseNoiseConfig, ...
%     AGCConfig=AGCConfig, ...
%     IqImbalanceConfig=IqImbalanceConfig);
%
% baseBandSignal = CPFSK(SampleRate=SampleRate, ...
%     TimeDuration=TimeDuration, ...
%     ModulatorOrder=ModulatorOrder, ...
%     SamplePerSymbol=SamplePerSymbol, ...
%     ModulatorConfig=ModulatorConfig);
%
% x1= baseBandSignal(x);
% x1 = txRF(x1);
% x1 = rayChannel(x1);
% y1 = rxRF({x1});
%
% %% 多天线场景
% NumTransmitAntennas = 1;
% NumReceiveAntennas = 2;
%
% txRF = TRFSimulator(StartTime=0.2, ...
%     SampleRate = SampleRate, ...
%     CarrierFrequency=CarrierFrequency, ...
%     IqImbalanceConfig=IqImbalanceConfig, ...
%     PhaseNoiseConfig=PhaseNoiseConfig, ...
%     MasterClockRate = MasterClockRate, ...
%     MemoryLessNonlinearityConfig=MemoryLessNonlinearityConfig);
% mimoChannel = MIMO(CarrierFrequency=CarrierFrequency, ...
%     SampleRate=SampleRate, PathDelays=0, ...
%     AveragePathGains=0, MaximumDopplerShift=0, ...
%     NumTransmitAntennas=NumTransmitAntennas, ...
%     NumReceiveAntennas=NumReceiveAntennas);
% rxRF = RRFSimulator(StartTime=0, TimeDuration=2, SampleRate=SampleRate, ...
%     NumReceiveAntennas=NumReceiveAntennas, CarrierFrequency=200e3, Bandwidth=20e3, ...
%     MasterClockRate=MasterClockRate, ...
%     MemoryLessNonlinearityConfig=MemoryLessNonlinearityConfig, ...
%     ThermalNoiseConfig=ThermalNoiseConfig, ...
%     PhaseNoiseConfig=PhaseNoiseConfig, ...
%     AGCConfig=AGCConfig, ...
%     IqImbalanceConfig=IqImbalanceConfig);
%
% baseBandSignal = CPFSK(SampleRate=SampleRate, ...
%     TimeDuration=TimeDuration, ...
%     ModulatorOrder=ModulatorOrder, ...
%     SamplePerSymbol=SamplePerSymbol, ...
%     NumTransmitAntennas = NumTransmitAntennas, ...
%     ModulatorConfig=ModulatorConfig);
%
% x2 = baseBandSignal(x);
% x2 = txRF(x2);
% x2 = mimoChannel(x2);
% y2 = rxRF({x2});

% Simple test for filter
M = 8; % Modulator order
cpfskMod = comm.CPFSKModulator(M, ...
    BitInput = false, SamplesPerSymbol = 32);
awgnChan = comm.AWGNChannel( ...
    NoiseMethod = 'Signal to noise ratio (SNR)', ...
    SNR = 0);
cpfskDemod = comm.CPFSKDemodulator(M, ...
    BitOutput = false, SamplesPerSymbol = 32);
spf = 10000; % Symobls per frame

delay = cpfskDemod.TracebackDepth;
errorRate = comm.ErrorRate( ...
    ReceiveDelay = delay);
data = randi([0 M - 1], spf, 1);
const = (- (M - 1):2:(M - 1))';
data = const(data(:) + 1);
modSignal = cpfskMod(data);

% Calculate occupied bandwidth
bw = obw(modSignal, SampleRate);
fprintf('Occupied Bandwidth: %.2f Hz\n', bw);

% Calculate theoretical bandwidth
h = 0.5; % modulation index
Rs = SampleRate / 32; % symbol rate (SamplesPerSymbol = 32)
theoretical_bw = (M - 1) * h * Rs + 2 * Rs; % Carson's rule for CPFSK
fprintf('Theoretical Bandwidth: %.2f Hz\n', theoretical_bw);

% Compare filtered signal
y = lowpass(modSignal, 14000, SampleRate, ImpulseResponse = "fir", Steepness = 0.99);
receivedData = cpfskDemod(modSignal);
errorStats = errorRate(data, receivedData);

receivedData1 = cpfskDemod(y);
errorStats1 = errorRate(data, receivedData1);
